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ABSTRACT 

We revisit the issue of spectral slope of MHD turbulence in the inertial range and 
argue that numerics favor Goldreich-Sridhar -5/3 slope rather than -3/2 slope. We 
also did precision measurements of anisotropy of MHD turbulence and determined 
the anisotropy constant Ca = 0.34 of Alfvenic turbulence. Together with previously 
measured Kolmogorov constant Ck = 4.2, or 3.3 for purely Alfvenic case, it constitutes 
a full description of MHD cascade in terms of spectral quantities, which is of high 
practical value for astrophysics. 

Key words: MHD - turbulence. 



1 INTRODUCTION 

Astrophysical plasmas are often described as ideal MHD 
fluid - a perfectly conducting, inviscid fluid described by 
MHD equations. The use of continuous fluid approach and 
absence of dissipation terms is motivated by the fact that 
astrophysical scales are typically humongous compared to 
the molecular and plasma mean free paths and the micro- 
scopic dissipation scalefl MHD can be successfully applied 
to all phases of the interstellar medium (ISM), large scales 
of molecular clouds, intracluster medium (ICM), stellar in- 
teriors and stellar outflows, etc. Reynolds numbers of as- 
trophysical flows are, typically, very high, which makes tur- 
bulence ubiquitous. Initially unmagnetized well-conductive 
turbulent fluid generates its own magnetic field which is dy- 
namically important on almost all scales. In spiral galaxies 
large-scale dynamo generates ordered fields on the scale of 
the galactic disc, while in more symmetric objects, such as 
galaxy clusters, small-scale dynamo operates and generates 
fields on scales up to 20 kpc. 

Inertial range of turbulence was introduced by 
lKolmogoro^l (|l94]J ) as a range of spatial scales where driv- 
ing and dissipation are unimportant and perturbations ex- 
ist due to energy transfer from one scale to another. In the 
inertial range of MHD turbulence perturbations of both ve- 
locity and magnetic field will be much smaller than the local 
Alfvenic velocity va = B/^/inp, due to the turbulence spec- 
trum being steeper than , therefore local mea n magnetic 
field will strongly affe ct dynamics in this range (|lroshnikovl 
ll964l : lKraichnanlll965h . As in the case of hydrodynamics, the 



^ This rule have some notable exceptions, such as molecular gas 
in the early Universe, which could be fairly viscous on scales of 
protogalaxies or the solar wind which is not collisional enough to 
fully support compressible modes. 



study of MHD turbulence began with weakly compressible 
and incompressible cases which are directly applicable to 
many environments, such as stellar interiors, ICM and hot 
phases of the ISM. Later it was realized that many features 
of incompressible MHD turbulence are still present even in 
supersoni c dynamics, due to the dominant effect of Alfv enic 
shearing (|Cho fc Lazaxiai i 2003: Bcre snvak et al.l 20051) . It 
had been pointed out by Goldreich &i Sridharl ( 19951 ) that 



strong mean field incompressible turbulence is split into the 
cascade of Alfvenic mode, described by Reduced MHD or 
RMHD (|Kadomtsev fc Pogutse|[l974l : IStrausslll976l ) and the 
passive cascade of slow (pseudo-Alfven) mode. In the strong 
mean field case it was sufficient to study only the Alfvenic 
dynamics, as it will determine all statistical properties of 
turbulence, such as spectrum or anisotropy. This decoupling 
was also observed in numerics. Luckily, being the limit of 
very strong mean field, RMHD has a precise two-parametric 
symmetry similar to the one in incompressible hydrodynam- 
ics (see §2), which, under certain conditions, makes universal 
cascade with power-law energy spectrum possible. 

Interaction of Alfvenic perturbations propagating in 
a strong mean field is unusual due to a peculiar disper- 
sion relation of Alfvenic mode, to = fc||DA, where fcy is a 
wavevector parallel to the mean magnetic field. This results 
in a tendency of MHD turbulence to create "perpendicu- 
lar cascade", where the flux of energy is preferentially di- 
rected perpendicular to the magnetic fleld. This tendency 
enhances the nonlinearity of the interaction, described by 
^ — 5vk_i_/vAk\i, which is the ratio of the mean-field term 
to the nonlinear term, and results in development of essen- 
tially strong turbulence. As turbulence becomes marginally 
strong, ^ ~ 1, the cascading timescales become close to the 
dynamical timescales Tease ~ f^yn = l/wk± and the per- 
turbation frequency uj has a lower bound due to an uncer- 
tainty relation Tcasci^ > 1 l|Goldreich fc Sridhaij[l995l ). This 
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makes turbulence being "stuck" in the ^ ~ 1 regime, which 
is known as "critical balance". Note, that the above argu- 
ment is not the only lower bound on uj, with another bound 
being due to the directional uncertainty of the va, which 
was discovered in iBeresnvak &: LazarianI (j2008;) . In this pa- 
per we will mostly consider so-called "balanced" turbulence 
in which both uncertainties coincide . We refer the reader to 
IBeresnvak fc LazarianI (|2008l . l2009bl ) and references therein 
for the more general imbalanced case. 

Goldreich-Sridhar model is predicting a k~^^^ energy 
spectru m with anisotrop y descr i bed as fcn ^ k^/^ . Numerica l 
studies iCho fc Vishniad (|2000l ): iMaron fc GoldreichI (120011 ) 
conf irmed steep spectrum and s cale-dependent anis o tropy , 
but iMaron fc GoldreichI l|200ll ): iMiiller fc GrappinI (|2005l ') 
claimed a shallower than —5/3 spectral slope in the strong 
mean field case, which was close to —3/2. T his motivated 
adjustments to the Goldreich-Sridhar m odel (jGaltier et al.l 
l2005l : lBoldvrev|[2005l : lGogoberidzj|2007l). A mode l with so 
called "dynamic aUgnment" ijBoldvrevlbOOSj . 120061 ) became 
popular after the scale-d ependent alignment was disc overed 
in numerical simulations (jBeresnvak fc Lazarian|[2006l ) . This 
model is based on the idea that the alignment between ve- 
locity and magnetic perturbations decreases the strength of 
the interaction scale-dependently, relying on the alignment 
being a power-law function of scale. This would, as they ar- 
gue, modify the spectral slope of MHD turbulence from the 
—5/3 Kolmogorov slope to the observed —3/2 slope. It was 
also claimed that there is a self-consistent turbulent mech- 
anism that produces such an alignment. In this paper we 
examine both the alignment and the spectrum. It turns out 
that the alignment is not universal and is tied to the driv- 
ing scale. Also, spectral resolution study at high numerical 
resolutions favor —5/3 spectral slope more than —3/2 slope. 
It seems that earlier measurements of the MHD slope were 
premature and were not conducted with enough rigor. 

In this paper we a) support our earlier claim of - 
5/3 scaling an d the measurement of Kolmogorov constant 
of [gcrcsnyak (l201ll ) and b) perform the measurement of 
anisotropy constant. These two measurements allow us to 
predict the local properties of MHD turbulence on any scale 
I, given only the dissipation rate e and Alfven velocity va- 
Such a straightforward prediction is of great practical use 
for astrophysics. 



2 BASIC EQUATIONS 

Ideal MHD equations describe the dynamics of ideally con- 
ducting inviscid fluid with magnetic field and can be written 
in Heaviside and c — 1 units as 



dtp + V-{pv) = 0, 
p(9t -I- V • V)v) = -VP + j X B, 
V • B = 0, 
9tB = V X (v X B), 

with current j = V x B and vorticity oj = V x v. This 
should be supplanted with energy equation and a prescrip- 
tion for pressure P. The incompressible limit assumes that 
the pressure is so high that the density is constant and veloc- 
ity is purely solenoidal (V ■ v = 0). This does not necessarily 



refer to the ratio of outer scale kinetic pressure to molecu- 
lar pressure, but could be interpreted as scale-dependent 
condition. Indeed, if we go to the frame of the fluid, local 
perturbations of velocity will diminish with scale and will 
be much smaller than the speed of sound. In this situation 
it will be possible to decompose velocity into low-amplitude 
sonic waves and essentially incompressible component of v, 
as long as we are not in the vicinity of a shock. The incom- 
pressible component, bound by V • v = 0, will be described 
by much simpler equations: 

dt-v = S{—uj X V 4- j X b), 
c»tb = V X (v X b), 

where we renormalized magnetic field to velocity units 
b = H/p^^^ (the absence of 47r is due to Heaviside units) 
and used solenoidal projection operator S — {1 — VA~^V) 
to get rid of pressure. Finally, in terms of Elsasser variables 

= V ± b this could be rewritten as 



dtw^ + ^(w^ ■ V)w^ 



0. 



(1) 



This equation resembles incompressible Euler's equa- 
tion. Indeed, hydrodynamics is just a limit of fe = in 
which w"*" — w~. This resemblance, however, is mislead- 
ing, as the local mean magnetic field could not be excluded 
by the choice of reference frame and, as we noted earlier, will 
strongly affect dynamics on all scales. We can explicitly in- 
troduce local mean field as va, assuming that it is constant, 
so that = w ± Va: 

dtSw^ =F (va ■ V)(5w* + SiSw"^ ■ V)(5w* = 0. (2) 

In the linear regime of small Sw^s they represent pertur- 
bations, propagating along and against the direction of the 
magnetic field, with nonlinear term describing their interac- 
tion. As we noted earlier, due to the resonance condition of 
Alfvenic perturbations they tend to create more perpendic- 
ular structure, making MHD turbulence progressively more 
anisotropic. This was empirically known from tokamak ex- 
periments and was used in so-called reduced MHD approx- 
imati on, which neglected parallel gradients in th e nonlinear 
term (jKadomtsev fc Poguts^ll974i : [Strausslll976l 'l . Indeed, if 
we denote || and _L as directions parallel and perpendicular 
to Va, the mean field term (waV||)(5u'* is much larger than 
{Sw^V \\)Sw'^ and the latter could be ignored in the inertial 
range where <C va- This will result in Equation [2] being 
split into 



dtSw., =F (va ■ V||)5w|| + ^((Jw]; • V±)5w|| 



dtSw^ =F (vA ■ V||)5wJ + S{Swf ■ Vx)5wj 



0, 



(3) 



(4) 



which, physically represent a limit of very strong mean 
field where 5wjj^ is a slow (pseudo-Alfven) mode and Sw"^ 
is the Alfven mode and Equation [3] describes a passive dy- 
namics of slow mode which is sheared by the Alfven mode, 
while Equation |4] describes essentially nonlinear dynamics 
of the Alfven mode and is known as reduced MHD. For our 
purposes, to figure out asymptotic behavior in the inertial 
range, it is sufficient to study Alfvenic dynamics and slow 
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mode can be always added later, because it will have the 
same statistics. 

It turns out that reduced MHD is often applicable be- 
yond incompressible MHD limit, in a highly coUisionless en- 
vironments, such as tokamaks or the solar wind. This is due 
to the fact that Alfven mode is transverse and does not re- 
quire pressure support. Indeed, Alfvenic perturbations rely 
on magnetic tension as a restoring force and it is sufficient 
that charged parti cles be tied to magn etic field lines to pro- 
vide inertia ( Schc kochihin et al.l|2009[ ) . 

A remarkable property of RMHD is that it has a precise 
two-parametric symmetry with respect to the anisotropy 
and the strength of the mean field; w — )■ wA, X — )■ XB, t — > 
tB/A, A — >■ KB/ A, which is similar to hydrodynamic sym- 
metry. Here A is a perpendicular scale, A is a parallel 
scale, A and B are arbitrary parameters of the transforma- 
tion. It is due to this precise symmetry and the absence 
of any designated scale, that we can hypoth esize univer- 
sal reg ime, similar to hydrodynamic cascade of iKolmogorovl 
(|l941 ). In nature, the universal regime for MHD can be 



achieved with <C «a. In numerical simulations, we can 
directly solve RMHD equations, which have precise symme- 
try already built in. From practical viewpoint, the statis- 
tics from the full MHD simulation with ~ O.Iwa is 
virtually indistinguishable from RMHD statistics and even 
5w^ ^ va is still fairly similar to the strong mean field case 
iBeresnvak fc Lazarianll2009bl ). 



3 BASIC SCALINGS 

As was shown in a rigorous perturbation study of weak 
MHD turbulence, it has a t endency of becoming stronger 
on smaller scales (jCaltier et al. 2000). Indeed, if fcy is con- 
stant and fcx is increasing, ^ = 5wk^/vAk\\ will increase, 
due to Sw ~ fc^ in this regime. This will naturally lead 
to strong turbulence, where ^ will stuck around unity due to 
two competing processes: 1) increasing interaction by per- 
pendicular cascade and 2) decrease of interaction due to the 
uncertainty relation TcascCJ > 1, where Tease is a cascad- 
ing timescale. Therefore, MHD turbulence will be always 
marginally strong in the inertial range, which means that 
cascading timescale is a ssociated with dynamical t imescale 
Tease ^ "^dyn — l/5'wki_ (|Goldreich fc Sridhaj [19951 ). In this 
case, assuming that energy transfer is local in scale and, 
therefore, depend only on perturbations amplitude on each 
scale, we can write Kolmogorov-type phenomenology as 



where is an energy fiux of each of the Elsasser vari- 
ables and is a characteristic perturbation amplitude on 
a scale A. Such an amplitude can be obtained by Fourier 
filterin g with a dyadic filter in k-space, see, e.g., ISeresnvaki 
l|2012l ). 

Since we consider so-called balanced case with both w's 
having the same statistical properties and energy fluxes, one 
of these equations is sufficient. This will result in a &w ~ 
A^''^, where A is a perpendicular scale, or, in terms of energy 
spectrum E{k), 



(6) 



where Ck is known as Kolmogorov constant. We will 
be interested in Kolmogorov constant for MHD turbulence. 
This scaling is supposed to work until dissipation effects kick 
in. In our further numerical argumentation dissipation scale 
will play a big role, but not from a physical, but rather from 
a formal point of view. We will introduce an idealized scalar 
dissipation term in a RHS of Eq. [T] as V^)"'^^w='=, 
where n is an order of viscosity and n = 2 correspond to 
normal Newtonian viscosity, while for n > 2 it is called hy- 
perviscosity. The dissipation scale for this Goldreich-Sridhar 
model is the same as the one for Kolmogorov model, i.e. 
rj — {u^/e)^^^^'^~^K This is a unique combination of Un and 
e that has units of length. Note that Reynolds number, es- 
timated as vL/v2, where L is an outer scale of turbulence, 
is around {L/t)Y^^ . 

Furthermore, the perturbations of w will be strongly 
anisotropic and this anisotropy can be calculated from the 
critical balance condition ^ ~ 1, so that fcy ~ k']^^ . Inter- 
estingly enough this could be obtained directly from units 
and the symmetry of RMHD equations from above. Indeed, 
in the RMHD limit, fcy or 1/A must be in a product with 
Va, since only the product enters the original RMHD equa- 
tions. We already assumed above that turbulence is local and 
each scale of turbulence has no knowledge of other scales, 
but only the local dissipation rate e. In this case the only 
dimensionally correct combination for the parallel scale A, 
corresponding to perpendicular scale A is 



A = CAVAX^'-'e 



2/3^-1/3 



(7) 



where we introduced a dimensionless "anisotropy con- 
stant" Ca- Equations [6] and [7] roughly describe the 
spectrum and anisotropy of MHD turbulence. Note, that 
Goldreich-Sridhar's —5/3 is a basic scaling that should be 
corrected for intermittency. This correction is negative due 
to structure function power-la w exponents being a concave 
function of their order (see, e.g. . iFrischI 1995f ) and is expected 
to be small in three-dimensional case. This correction for 
hydrodynamic turbulence is around —0.03. Such a small de- 
viation should be irrelevant in the context of debate between 

—5/3 and —3/2, which differ by about 0.17^ 

An alternative model was proposed by iBoldvrevI (|2005l . 

who suggested that these scalings are modified by a 
scale dependent factor that decreases the strength of the in- 
teraction, so that RHS of the Equation [5] is effectively mul- 
tiplied by a factor of {I/LY^^ , where L is an outer scale. We 
will discuss this hypothesis further in §8. In this case the 
spectrum will be expressed as E{k) = CK2(?^^k~^^'^ L^'^ . 
Note that this spectrum is the only dimensionally correct 
spectrum with k~^^^ scaling, which does not contain dissi- 
pation scale rj. The absence of L/rj, is due to so-called zeroth 
law of turbulence which states that the amplitude at the 
outer scale should not depend on the viscosity. This law fol- 
lows from the locality of energy transfer has been know em- 
pirically to hold very well. The dissipation scale of Boldyrev 
model is different from that of the Goldreich-Sridhar model 
and can be expressed as r?' = (t,3/^)i/{3n-i.5)^o.5/(3n-i.5) ^ 
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Table 1. Three-dimensional RMHD simulations 



Run Tlx ■ riy ■ Dissipation (e) L/rj 

256 ■ 768^ -6.82 • lO-^^fc** 0.073 200" 

R2 512 • 1536^ -1.51 • lO-i^feS 0.073 400 

R3 1024 ■ 3072^ -3.33 ■ 10~^'^fc'=' 0.073 800 

'RA 76p -6.82 ■ lO^^fc** 0.073 200" 

R5 ISSO'^ -1.51 ■ lO-^^fc'^ 0.073 400 

~Rj6 384 • 1024^ -1.70 ■ 10-*fc^ 0.081 280" 

R7 768 ■ 2048^ -6.73 ■ IQ-^k^ 0.081 560 

^!8 76p -1.26 ■ 10-*A;^ 0.073 350" 

R9 1536=* -5.00 ■ 10^5^2 o.073 700 



been obtained with rather moderate resolutions, which im- 
plies that hydrodynamic cascade has fairly narrow locality. 
When the resolution is increasing, the convergence is sup- 
posed to become better. This is because of the increased sep- 
aration of scale between driving and dissipation. Better con- 
vergence means higher precision in d iscriminating differen t 
scalings, which was demonstrated in iKaneda et all l|2003h . 
which directly measured the intermittency correction to the 
spectral slope. The optimum strategy for our MHD study 
is to perform the largest resolution simulations possible and 
do the scaling study as described in this section. 



4 THE SCALING ARGUMENT 

Turbulence with very long ranges of scales is common in 
astrophysics. Three-dimentional numerics, however, are un- 
able to reproduce such ranges and normally struggles to ob- 
tain even a small inertial range. In this situation we have to 
use rigorous quantitative arguments in order to investigate 
asymptotic scalings. 

Suppose we performed several simulations with different 
Reynolds numbers. If we believe that turbulence is univer- 
sal and the scale separation between forcing scale and dis- 
sipation scale is large enough, the properties of small scales 
should not depend on how turbulence was driven. This is 
because neither MHD nor hydrodynamic equations explic- 
itly contain any scale, so simulation with a smaller dissipa- 
tion scale could be considered, due to symmetries from §2, 
as a simulation with the same dissipation scale, but larger 
driving scale. E.g., the small scale statistics in a 1024^ simu- 
lation should look the same as small-scale statistics in 512^, 
if the physical size of the elementary cell is the same and 
the dissipation scale is the same. 

This scaling argument requires that the geometry of 
the elementary cells is the same and the actual numerical 
scheme used to solve the equations is the same. Also, nu- 
merical equations should not contain any scale explicitly, 
which is usually satisfied. The scaling argument does not re- 
quire high precision on the dissipation scale or a particular 
form of dissipation, either explicit or numerical. This is be- 
cause we only need the small-scale statistics to be similar in 
both simulations. This is achieved, on one hand, by applying 
the same numerical procedure, and, on the other hand, by 
turbulence locality that ensures that outer scale influence is 
small. 

In practice, the scaling argument also known as a res- 
olution study is done in the following way: the averaged 
spectra in two simulations are expressed in dimensionless 
units corresponding to the expected scaling, for example a 
E{k)k^^^e~'^^^ is used for hydrodynamics, and plotted ver- 
sus dimensionless wavenumber kr^, where dissipation scale r; 
correspond to the same model, e.g. t] = {v^/e)^^'^ is used 
for scalar second order viscosity f and Kolmogorov phe- 
nomenology. As long as the spectra are plotted this way and 
the scaling is correct, the curves obtained in simulations with 
different resolutions has to converge on small scales. Not 
only the spectrum, but any other statistical property of tur- 
bulence can be treated this way. This metho d has been used 
in hyd r odynamics s ince long time ago, e.g. bv'Y eung fc Zhoul 
l|l997h : iGotoh et al. (2002); Kaneda et al.. (,200^ ') For hv- 
drodynamics good convergence on the dissipation scale has 



5 NUMERICAL METHODS 

We used pseudospectral dealiased code to solve RMHD 
equations. Same code was used earlier for RMHD, incom- 
pressible MHD and incompressible hydrodynamic simula- 
tions. The RHS of Eq. |4]was complemented by an explicit 
dissipation term — V^)"''^w* and forcing term f. The 
code and the choice for numerical resolution, driving, etc, 
was desc ribed in g reat detail in our earlier p ublications 
( Bcrcsnv ak fc Lazar ian 2009bllaL [2OI0I : iBeresnvak 2011). Ta- 
ble 1 shows the parameters of the simulations. The Kol- 
mogorov scale is deflned as 7; = (i/j^/e)^'''-^""^' , the integral 
scale L — iv/AE JJ" k^^E{k) dk (which was approximately 
0.79 for Rl-3). Dimensionless ratio L/rj could serve as a 
"length of the spectrum" , although spectrum is actually sig- 
nificantly shorter for n=2 viscosity and somewhat shorter for 
n=6 viscosity. 

Since we would like to use this paper to illustrate the 
resolution study argument we used a variety of resolution, 
dissipation and driving schemes. There are four schemes, 
presented in Table 1, and used in simulations Rl-3, R4-5, 
R6-7 and R8-9. In some of the simulations the resolution in 
the direction parallel to the mean magnetic field, Ux, was 
reduced by a factor compared to perpendicular resolution. 
This was deemed possible due to an empirically known lack 
of energy in the parall el direction in fc-space an d has been 
used before (see, e.g., iMiiUer fc Grappinll200d ). The R4-5 
and R8-9 groups of simulations were fully resolved in paral- 
lel direction. One would expect that roughly the same res- 
oluti on will be required in parallel and perpendicular direc- 
tion (jBeresnvak fc Lazarian|[2009bl ) . In all simulation groups 
time step was strictly inversely proportional to the resolu- 
tion, so that we can utilize the scaling argument. 

Driving had a constant energy injection rate for all sim- 
ulations except R6-7, which had fully stochastic driving. All 
simulations except R8-9 had Elsasser driving, while R8-9 
had velocity driving. All simulations were well-resolved and 
R6-7 were overresolved by a factor of 1.6 in scale (a fac- 
tor of 2 in Re). The anisotropy of driving was that of a 
box, while injection rate was chosen so that the amplitude 
was around unity on outer scale, this roughly corresponds 
to critical balance on outer scale. Indeed, as we will show 
in subsequent section, since anisotropy constant is smaller 
than unity, our driving with A ~ A ~ 1 and Sw ~ 1 on outer 
scale is somewhat over-critical, so A decreases after driving 
scale to satisfy uncertainty relation (see Fig. [5]). This is good 
for maintaining critical balance over wide range of scales as 
it eliminates possibility for weak turbulence. 

In presenting four groups of simulations, with different 
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Figure 1. Numerical convergence of spectra in all simulations. Two upper rows are used to study convergence assuming Boldyrev 
model and two bottom rows - assuming Goldreich-Sridhar model. Note that definition of dissipation scale r] depends on the model, see 
§3, this difference is tiny in hyperviscous simulations Rl-5, but significant in viscous s imulations R 6 -9. N umerical convergence require 
that spectra will be similar on small scales, including the dissipation scale, see, e.g. iGotoh et al As we see from the plots, 

numerical convergence is absent for Boldyrev model. For Goldreich-Sridhar model the convergence is reached only at the dissipation 
scale. Higher-resolution simulations are required to demonstrate convergence in the inertial range. 



geometries of elementary cell, different dissipation terms and 
different driving, our intention is to show that the scaling 
argument works irrespective of numerical effects, but rather 
relies on scale separation and the assumption of universal 
sc aling. Simu l ations Rl-3 are the same as those presented 
in iBeresnvakI (|201l h. Largest simulation R3 with 3.1 x 10^ 
points used about 0.85 million CPU hours on TACC Ranger. 



6 DRIVING AND CONSERVATION LAWS 

All of our simulations except R8-9 used Elsasser driving. 
The choice of this driving is due to the formulation of our 
problem - we want to obtain asymptotic scalings in the in- 
ertial range of MHD turbulence. This is also the reason we 
use RMHD, although sometimes full MHD with strong mean 
fi eld is used ( MuUer fc Grappin 2005; Bcrcsnyak & Lazarian 
[2009ai lbl). In this situation we simulate a tiny, compared to 
the outer scale, volume of turbulence. Since incompressible 
MHD has two Elsasser energy cascades and exchange of en- 
ergy between these cascades is impossible, it make sense that 
we provide energy to these two cascades independently, in 
fact this is the onl y way to simulate imbalanced turbulence 
in a periodic box (jBeresnvak fc Lazaria nl l2009bl ). 

When driving is concerned, often the issue of integral 
conservation laws is raised. In particular, inviscid hydro- 
dynamics has Kelvin's circulation theorem which preserves 
circulation along any closed curve, moving with the fluid. 



as long as the external force is potential. This also implies 
that the lines of vorticity are "frozen" into the fluid. Natu- 
rally, driving force in most hydrodynamic simulations is not 
potential and, as a result, Kelvin's theorem is broken. Al- 
though non-potential forces are known to drive turbulence 
in nature, this usually is only important on the outer scale 
of turbulence, but in the inertial range such forces should be 
negligible compared to the dynamic pressure. 

What is the physical meaning of volumetric force that is 
typically used in simulation s which try to reproduce inertial 
range, such as lGotoh et al.l l|2002l )? The key to understand 
this is to remember that each fluid element is embedded into 
the surrounding fluid that exert forces on its boundaries. 
These forces supply energy that drives turbulence within 
this element. Similarly, in MHD fluid an external electro- 
motive force (EMF) exists on the boundaries of the fluid 
element and the curl of this EMF enters in RHS of the in- 
duction equation. Unfortunately, simulating fluid element 
embedded into a larger system would effectively mean simu- 
lating this larger system and is not practical. In practice we 
use periodic boxes to mimic statistical homogeneity. Despite 
driving in such boxes being restricted to low wavenumbers, 
the correlation scale is typically smaller than the box size 
by a factor of 3-4, which simulates many fluid elements in 
a box. The driving around k=2..3, therefore, emulates an 
action exerted on the boundaries of these fluid elements. If 
fcmax is a maximum wavenumber of driving, then the change 
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Figure 2. Resolution study for "dynamic alignment", assuming 
Boldyrev scaling. Both axis are dimensionless, solid is higher res- 
olution and dashed is lower resolution. Convergence is absent for 
all simulations. This suggests that l^ '^^ is not a universal scaling 
for alignment. 



of circulation within a loop of size I < 27r/femax, associated 
with forcing, is going to be small, emulating the situation of 
real turbulence where the forces exerted on the boundaries 
of fluid element will actually preserve circulation along any 
loop within this element. 

Similarly, Alfven's theorem is a conservation of the mag- 
netic flux through a given loop and this implies that the 
magnetic field lines are frozen into the fluid. It can be bro- 
ken by an external EMF, however, if such an EMF is re- 
stricted to low wavenumbers in a simulation with periodic 
box, this will effectively emulate an external EMF action on 
the boundaries of the fluid element. There is a full analogy 
between Alfven's theorem and Kelvin's theorem since we 
expect neither significant external forces nor external EMF 
at the inertial range scales. Also, it turns out that there is 
an analogy between the breakdown of Kelvin's theorem and 
Alfven's theorem in turbulent flows (|Evink fc Aluiell2006l : 
lEvinkll2007l ). 

To summarize, the driving for a simulation with a strong 
mean field or RMHD simulation reproducing inertial range, 
must have Elsasser driving, i.e. independent driving of w"*" 
and . This will simulate supply of Elsasser energies from 
larger scales. If turbulence is strong, a pure velocity driv- 
ing is also possible due to the quick nonlinear decorrelation 
of and w~ eddies. Volumetric Elsasser driving does not 
honor Kelvin's or Alfven's theorems, however, by restricting 
driving to large scales we effectively emulate the action of 
external forces which conserve fluxes and circulations. On 
the dissipation scales Kelvin's or Alfven's theorems are bro- 
ken by viscosity and magnetic diffusivity and in the inertial 
range they will be broken by turbulence. Therefore, there 
is an analogy between forced viscous hydrodynamic simula- 
tions in a periodic box and forced dissipative MHD simula- 
tions. 



7 SPECTRA 

Fig. [1] presents a resolution study all simulations. The upper 
rows assume Boldyrev scaling, while the bottom rows as- 
sume Goldreich-Sridhar scaling. Reasonable convergence on 
small scales was achieved only for Goldreich-Sridhar scaling. 
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Figure 3. DA slope, defined as l/DAdT)A/dl, solid is higher res- 
olution and dashed is lower resolution. Dynamic alignment slope 
does not converge and has a tendency of becoming smaller in 
higher-resolution simulations. This may indicate that the asymp- 
totic alignment slope is zero, which will correspond to the GS95 
model. 



The normalized amplitude at the dissipation scale for two 
upper rows of plots systematically goes down with resolu- 
tion, suggesting that —3/2 is not an asymptotic scaling. The 
flat part of the normalized spectrum on Rl-3 plots was fit 
to obtain Kolmo gorov const a nt oi Cka = 3.27 ± 0.07 which 
was reported in iBeresnvakI |20l3). The total Kolmogorov 
constant for both Alfven and slow mode in the above paper 
was estimated as Ck = 4.2 ±0.2 for the case of isotropically 
driven turbulence with zero mean field, where the energy 
ratio of slow and Alfven mode C's is between 1 and 1.3. This 
larger value Ck = Ca'a(1 + Ca)^^^ is due to slow mode be- 
ing passively advected and not contributing to nonlinearity. 
The measurement of Cka had relied on an assumption that 
the region around krj ~ 0.07 represent asymptotic regime. 
It is possible, though unlikely, that Cka is slightly under- 
estimated due to this region is somewhat lower than the 
asymptotic regime due to antibottleneck effect exacerbated 
by a reduced parallel resolution. The antibottleneck effect, 
however is typically much smaller that the bottleneck effect, 
so this correction is probably within the stated errorbars. 
Next simulation in a R4-5 series (3072=*) should make this 
clear. As far as simulations with normal viscosity R6-9 go, 
it seems impossible to approach good inertial range until 
resolutions of at least 4096"^. 



8 DYNAMIC ALIGNMENT 



iBoldvrevI (|2005l ) suggested that and w~ eddies are sys- 
te matically aligned. This was investigated numerically in 
in iBeresnvak fc LazarianI (|2006l ) and no significant align- 
ment was found for the averaged angle between w"*" and 
w~, AA = {|5w^ X 5w^|/|5wjJ'||5w^|), but when this 



angle was weighted with the amplitude PI 



1) /(|( 5 w^ I |(5w 7 1), some alignment was found. Then 
iBoldvrevI (|2006|) sugge sted alignment between v and b 
and ( Mason et al.ll2006l ) suggested a particular amplitude- 
weighted measure, DA = {\5vx x 5bA|)/{|5vA||5bA|). We 
note that DA is similar to PI but contain two effects: 
alignment and local imbalance. The latter could be mea- 
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sured with IM = mwt? - S(w- f\)/(5(w+f + 5{w-f), 

iBeresnvak fc LazariarJ (|2009al ) . 

In this section we check the assertion of lBoldvrevI l|2005l . 

I2OO6I ) that alignment depends on scale as A^'^*, by using DA 
which is, by some reason, favored by aforementioned group. 
We did a resolution study of DA, assuming suggested scal- 
ing, which is presented on Figure [2] Convergence was absent 
in all simulations. Note, that previous studies that claimed 
that there is a good correspondence with Boldyrev model 
did not perform the scaling study, therefore these claims 
are not well substantiated. A result from a single isolated 
simulation could be easily contaminated by the effects of 
outer scale, since it is not known a-priori how local MIfD 
turbulence is and what resolution is sufficient to get rid of 
such effects. On the contrary, the resolution study offers a 
systematic approach to this problem. 

Fig. [3] shows "dynamic alignment " slope for all simula- 
tions. Although there is no convergence as in the previous 
plot, it is interesting to note that alignment slope decreases 
with resolution. This suggests that most likely the asymp- 
totic state for the alignment slope is zero, i.e. alignment 
is scale-independent and Goldreich-Sridhar model is recov- 
ered. Also, alignment from simulations Rl-5 seems to indi- 
cate that the maximum of the alignment slope is tied to the 

outer scale, therefore alignme nt is a transitional effec t. 

In our earUer studies (|Beresnvak fc LazariarJ |2006| . 

l2009ah we measured several types of alignment and 
found no evidence that all alignment measures follow the 
same scaling, see, e.g., Fig. [H As one alignment mea- 
sure, PI, has been already known to be scale-dependent 
IBeresnvak fc Lazarian|[2006l ) prior to DA, it appears that a 
particular measure of the alignment in iMasoii et all (|2006l ') 
was been picked for being most scale-dependent and no thor- 
ough explanation was given why it was preferred. Futher- 
more, it was claimed that DA has the asymptotic scaling of 
~ (Z/L)^'''* and at the same time it is an interaction weak- 
ening factor that will result in a —3/2 spectrum. This is 
unlikely, since the interaction weakening factor will appear 
from a third-order structure function, while DA is a correc- 
tion factor based on a ratio of second-order structure func- 
tions. Intermittency corrections are supposed to be small, 
however. More importantly, a physical justification of DA 
and its preference over other measures, such as PI or IM, 
that differ from DA quite significantly, was lacking. 

We are not aware of any convincing physical argumen- 
tation explaining why alignment s hould necessarily be a 
power-law of scale. iBoldvrevI l|2006l ) argues that alignment 
will tend to increase, but will be bounded by field wander- 
ing, i.e. the alignment on each scale will be created indepen- 
dently of other scales and will be proportional to the rel- 
ative perturbation amplitude SB/B. But this violates two- 
parametric symmetry of RMHD equations mentioned above, 
which suggests that field wandering can not destroy align- 
ment or imbalance. Indeed, a perfectly aligned state, e.g., 
with 5w~ = is a precise solution of MHD equations and it 
is not destroyed by its own field wandering. The alignment 
measured in simulations of strong MHD turbulence with dif- 
ferent values of SB l/Bq showed very little or no d ependence 
on this parameter l|Beresnvak fc Lazaria nl l2009al ). 

Why alignment measures are scale-dependent quanti- 
ties over about one order of magnitude in scale is an inter- 
esting question. Most plausible explanation is that because 
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Figure 4. Slopes of several alignment measures vs scale in R4-5 
(for definitions see the text). Each measure follows its own scaling, 
however there are indications that they are all tied to the outer 
scale, due to the maximum of alignment being a fraction of the 
outer scale, which is an indication that their scale-dependency is 
of transient nature. 



MHD turbulence is much less local t han hydro turbu lence 
(|Beresnvak fc LazarianI l2009al . l2010l : iBeresnvakI HoU) and 
because driving does not mimic the properties of the iner- 
tial range, the transition to asymptotic statistics is very wide 
and many quantities appear as scale-dependent, while they 
are simply adjusting to the asymptotic regime. 

The contribution to energy flux from different k wave- 
bands is important to understand, since most cascade mod- 
els assume locality, or rather to say the very term " cascade" 
assumes locality. An analytical upper bound on locality sug- 
gests tha t the width of th e energy transfer window can scale 
as l|Beresnvakl |2012| ) . however, in practice turbulence 

can be more local. The observation of lBeresnvak fc LazarianI 
(j2009ah that MHD simulations normally lack bottleneck ef- 
fect, even with high-order dissipation, while hydrodynamic 
simulations always have bottleneck, which is especially dra- 
matic with high-order dissipation, is consistent with above 
conjecture on locality, since bottleneck effect relies on local- 
ity of energy transfer. As locality constraint depends on the 
efficiency of the energy transfer, so that the efficient energy 
transf er must be local, w hile inefficient one could be no nlo- 
cal (Ber esnvak fc Lazaria n 2010; Bcrcsnyak 20Ill l2012h . As 
we observe larger Ck in MHD turbulence compared to hy- 
drodynamic turbulence, the former could be less local than 
the latter, which is consistent with our earlier findings. 



9 ANISOTROPY 

In Section 3 we suggested that anisotropy should be 
universal in the inertial range and expressed as A = 
CA«AA^/3e-'/^ where Ca is an anisotropy constant to 
be determined from the numerical experiment or obser- 
vation. Note, that both Alfvenic and slow modes should 
have the same anisotropy. This is because they have the 
same ratio of propagation to nonlinear timescales. Fig- 
ure [5] shows anisotropy for the two best resolved groups 
Rl-3 and R4-5. We used a model independent method of 
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Figure 5. The scaling study for anisotropy shows moderately 
good convergence to a universal anisotropy A = Ca'''a^^^^^~^^^ 
with anisotropy constant Ca of around 0.34. 



minimum parallel structure fu nction, described in detail in 
iBeresnvak fc LazarianI (|2009lJ ) . Alternative definitions of lo- 
cal mean field give comparable results, as long as they are 
reasonable. The convergence of anisotropy curves at the dis- 
sipation scale is not as good as convergence of spectra in 
Rl-5 on the dissipation scale. This is due to this measure 
being calculated from structure functions which are less local 
in scale than 3D spectra, i.e. disipation scale is still being 
somewhat affected by the transition to asympotic regime. 
From Rl-3 we obtain Ca = 0.34. Note, that the conven- 
tional definition of critical balance involve the amplitude, 
rather than (eX)^''^, so the constant in this classical for- 

1/2 

mulation will be CaCj^ « 0.63, which is closer to unity. 
Together with energy spectrum this is a full description of 
universal axisymmetric two-dimensional spectrum of MHD 
turbulence in the inertial range. 



10 SUMMARY AND DISCUSSION 

In this paper we argued that the properties of Alfven and 
slow components of MHD turbulence in the inertial range 
will be determined only by the Alfven speed va, dissipation 
rate e and the scale of interest A. The energy spectrum and 
anisotropy of Alfven mode will be expressed as 



E{k) = CKe 



2/3 I 



5/3 



A/\ = CAVA{\e)-'^\ 

with Ck ~ 3.3 and Ca = 0.34. If the slow mode is 
present, its anisotropy will be the same, and it will con- 
tribute to both energy and dissipation rate. Assuming the 
ratio of slow to Alfven energies between 1 and 1.3, the latter 
was observed in statistically isotropic high resolution MHD 
simulation with no me an field, we can u se Ck ~ 4.2 for the 
total energy spectrum (|Beresnvakll201ll '). 

Anisotropy of MHD turbulence is an important prop- 
erty that affects such processes as interaction with cosmic 



rays, see, e.g., lYan fc LazarianI (|2002l) . bmce cosmic ray pres- 
sure in our Galaxy is of the same order as dynamic pressure, 
their importance should not be underestimated. Another 
process affected is the three-dimentional tu rbulent recon- 
nection, see, e.g., iLazarian fc Vishnia3 |l99^. 

Previous measurements of the energy slope relied on the 
highest-resolution simulation and fitted the slope in the fixed 
fc-range close to the driving scale, typically between A; = 5 
and k = 20. We argue that such a fit is unphysical unless a 
numerical convergence has been demonstrated. We can plot 
the spectrum vs dimensionless k-q and if we clearly see a 
converged dissipation range and a bottleneck range, we can 
assume that larger scales, in terms of krj represent inertial 
range. In fitting fixed fc-range at low k we will never get 
rid of the influence of the driving scale. In fitting a fixed krj 
range, the effects of the driving will diminish with increasing 
resolution. 

Since we still have trouble transitioning into the inertial 
range in large mean field simulations, for now it is impos- 
sible to demonstrate inertial range in st atistically isotropic 
simula tions similar to once presented in lMiiller fc Grappinl 
(|2005l ). This is because we do not expect a universal power- 
law scaling in transAlfvenic regime, due to the absence of 
appropriate symmetries and the transitioning to subAlfvenic 
regime, where such scaling is possible, will require some ex- 
tra scale separation. These two transitions require numeri- 
cal resolution that is even higher than the highest resolution 
presented in this paper and for now seem computationally 
impossible. 

Full compressible MHD equations contain extra degrees 
of freedom, which, in a weakly compressible case, entails the 
additional cascade of the fast MHD mode, possibly of weak 
natur e. Supersonic simulati ons with moderate Mach num- 
bers (|Cho fc LazarianI l2003h show that Alfvenic cascade is 
pretty resilient and is not much affected by compressible 
motions. The models of the "universal" supersonic turbu- 
lence covering supersonic large scales and effectively sub- 
sonic small scales are based mainly on simulations with lim- 
ited resolution and unlikely to hold true. This is further rein- 
forced by the results of this paper which demonstrated that 
even a much simpler case of sub- Alfvenic turbulence require 
fairly high resolutions to obtain an asymptotic scaling. 

In this paper we treated so called balanced case with 
(Jul"*" ~ 5w ~ . A more general imbalanced case has been 
discussed in IBeresnvak fc LazarianI (|2008l . l2009bl . l2010f) . see 
also references therein. 



10.1 Results of iGrappin &: Miilleil (|2010l ^ 

iGrappin fc MiiUeil (|2010l ) observed that MHD turbulence 
have scale-independent anisotropy, in contrast with the 
Goldreich-Sridhar model and our own results. We believe 
that this measurement is correct, in fact, similar 2D spec- 
tra has been reported in our study IBeresnvak fc LazarianI 
(i2009bl ). However, these are 2D spectra obtained with re- 
spect to the global mean magnetic field . A triv ial exercise 
ICho et all |2002| : IBeresnvak fc LazarianI 1200931 ) show that 
measuring anisotropy with respect to the global field will 
destroy scale-dependent anisotropy, even in the case of very 
strong field. Indeed, even if we have SBl/Bq ^ 1, the 
anisotropy in the global frame will be limited from above by 
Bq/SBl, due to field wandering, while the Goldreich-Sridhar 
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anisotropy on the scale / will be much higher, ~ Bq/SBi, by 
a factor of Bl/Bi. In particular, each individual volume on 
scale I will have anisotropy with respect to the mean field 
averaged on the same volume, however, the direction of this 
field will deviate from the global mean field statistically, with 
RMS deviation of around SBl/Bo- This directional devia- 
tions will easily destroy anisotropy ~ Bq/SBi if we aver- 
age spectra measured with respect to Bo- Therefore, one 
must measure anisot ropy with respec t to lo cal mean field, 
which was realized in lCho fc Vishniad (|2000l '). From this ar- 
gument, we see that the local mean field must be averaged 
on scales smaller or equal to I, a scale at which spectrum 
or structure function is measured and I itself is a preferred 
averaging scale. Also, it is anisotropy with respect to the lo- 
c al mean field, which is imp ortant for cosmic ray dynamics 
(Yan fc Lazarianll2002l . l2004l ). 



10.2 Results of lMason et al.l (|201lh 

iMason et alj (|201ll ) conducted a series of low-resolution sim- 
ulations (256^ — 512^) studying the effect on dynamic align- 
ment by numerics on viscous scale, e.g., by changing resolu- 
tion with constant Re, or changing elementary cell geome- 
try, such as the ratio of parallel to perpendicular resolution. 
Their conclusion was that alignment is, indeed, somewhat 
infiuenced by numerics. 

As we explained previously, numerical effects on vis- 
cous scale are irrelevant for scaling study, if latter was 
conducted p roperly. Therefore it was quite puzzling that 



Mason et al.l (20 11') claimed that their results invalidate 



Beresny^iioTJ), which did a proper scaling study with 



high resolution simulatio ns. Furthermore, it is also puzzling 
that ason et al] (|201ll ) claimed an "extended scaling law" 
for alignment (previously mentioned A^^*), despite the fact 
that they did not even attempt a proper scaling study for 
this quantity and did not demonstrate convergence. 

Although studying obscure effects on viscous scales is 
pretty common in modern numerical studies of hydrody- 
namic turbulence, we will argue that such studies in MHD 
are of limited value for astrophysics, since dissipation mech- 
anism of astrophysical plasma is, typically, not a scalar 
second-order diffusion. 

On the contrary, an asymptotic scalings of the iner- 
tial range of MHD turbulence is of high value. However, 
such scalings could only be demons trated by a proper rig - 
orous scaling study in the spirit of I Yeung fc Zhoul (|l997h : 
iGotoh et al.r(|2002h . and not by arbitrarily assigning "iner- 
tial range" to a fixed range in wavenumbers. 
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